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ABSTRACT 


Jentink, Thomas Neil. MSAA, Purdue University, August, 1989. Formulation Of 
Boundary Conditions For The Multigrid Acceleration Of The Euler And Navier Stokes 
Equations. Major Professor: William J. Usab, Jr. 

An explicit, Multigrid algorithm has been written to solve the Euler and Navier 
Stokes equations with special consideration given to the coarse mesh boundary 
conditions. These are formulated in a manner consistent with the interior solution, 
utilizing forcing terms to prevent coarse-mesh truncation error from affecting the fine- 
mesh solution. A 4-Stage Hybrid Runge-Kutta Scheme is used to advance the solution 
in time, and Multigrid convergence is further enhanced by using local time-stepping 
and implicit residual smoothing. Details of the algorithm are presented along with a 
description of Jameson’s standard Multigrid method and a new approach to formulating 
the Multigrid equations. This approach utilizes a general filtering operator to derive the 
coarse-mesh equations in partial differential form, and is reduced to Jameson’s 
Multigrid scheme by specifying a particular discrete filter. The correct boundary 
transfer operator is formulated from this filter, and forcing terms are then derived for 
the coarse-mesh boundary conditions. Results are presented for inviscid channel flow 
of subsonic, transonic, and supersonic speeds over a circular-arc bump, viscous flat plat 
flow, viscous flow over a circular-arc bump, and inviscid and viscous flow over a VKI 
gas turbine rotor blade. These results will show the importance of the correct 



implementation of the coarse-mesh boundary conditions by comparison of convergence 
levels with and without the boundary forcing terms that were derived. 



1 


CHAPTER 1 • INTRODUCTION 


Advances in computer technology over the past few years , coupled with the rising 
cost of experimentation has resulted in the increased role of computational methods in 
the design process. The improved algorithm development occurring along with these 
advances has made it possible to surpass the level of approximation provided by 
potential solvers. The higher order Euler and Navier-Stokes equations can now be 
solved for a number of more realistic problems covering a wide range of 
complexity. t 1-4 ' Finite-volume, explicit time-stepping schemes have been used quite 
extensively in solving these problems for both the Euler and Navier-Stokes 
equations^ 1-9 ' Explicit schemes are generally favored over implicit schemes due to 
their ease of implementation, minimum number of computations, and low storage 
requirements. Although explicit schemes have been successful, their main disadvantage 
over implicit schemes is their stability limit. Whereas an implicit scheme is by theory 
unconditionally stable and therefore able to use very large time-steps, an explicit 
scheme can be quite limited by the maximum allowable time-step that may be used. In 
order to retain the advantages of an explicit scheme, while still reaching solutions 
efficiently, a large amount of research has focused on improving the performance of 
explicit time-stepping schemes. 
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As part of this effort, a class of hybrid Runge-Kutta time-stepping schemes was 
developed by Jameson.^ Runge-Kutta schemes are usually chosen for their high level 
of time accuracy, but in the interest of steady-state solutions, Jameson neglected time 
accuracy and tuned the Runge-Kutta coefficients to provide the largest possible time 
step and the best error damping characteristics. To further enhance the convergence of 
steady-state solutions, three standard acceleration methods are also used: 

— Local time-stepping 

— Implicit residual smoothing 

— Multiple-Grid and Multigrid Acceleration 

The time-step for a given computational cell is proportional to its size. 
Consequently, the convergence of a scheme may be severly limited by mesh geometry 
rather than the actual physics of the problem if very small mesh spacings are present, 
since the smallest time step is the limiting factor in global time-stepping (required for 
time-accurate problems). For steady-state solutions, since time-accuracy is not 
required, the solution in each cell may be marched in time by its local time-step, greatly 
increasing convergence rates. 

Implicit residual smoothing is a method introduced by Jameson^ 61 for the Euler 
equations in which the stability limit (CFL number) may be increased to a large degree 
by replacing the residual at each point by a implicit average based on the residuals in 
the rest of the domain. This smoothing, used in conjunction with Jameson’s hybrid 
Runge-Kutta schemes can increase the maximum allowable CFL number by three times 
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over that given by the stability limit. ^ 

Multigrid is an acceleration tool first developed by Achi Brandt^ 10, ^ to increase 
convergence rates for elliptic-type problems. In elliptic problems, high frequency 
errors are eliminated quickly by the relaxation scheme. The limiting factor on 
convergence is the low frequency errors that remain. In order to more rapidly eliminate 
these low frequency errors, Brandt utilized a sequence of meshes made up of 
successively larger cells. Relaxation sweeps are made on the initial, finest mesh, 
effectively damping the high frequency errors (wavelengths on the order of the mesh 
spacings). Then, the solution is transferred by an appropriate interpolation scheme, 
called a transfer operator, to the next coarser mesh, where the relaxation sweeps are 
again performed. Now, since the cells are larger, the relaxation scheme removes error 
frequencies that are lower (longer wavelengths) with respect to the fine mesh solution. 
When the solution is transferred back to the fine mesh, a large amount of the low 
frequency errors have been eliminated by the relaxation on the coarse mesh levels. On 
the coarse levels, forcing terms (source terms added to the coarse-mesh equations) are 
utilized to prevent the coarse-mesh accuracy from affecting the level of accuracy of the 
final fine-mesh solution. Brandt’s Multigrid scheme is labelled the Full Approximation 
Storage (FAS) scheme due to the fact that the full solution from the finest mesh is 
transferred to the coarser mesh levels. 

The first application of this type of acceleration scheme for the hyperbolic Euler 
equations was introduced by Ni . [121 His Multiple-Grid acceleration scheme differs 
from Brandt’s FAS Multigrid scheme in that only changes in the fine-mesh solution. 
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versus the full solution, are transferred to the coarse mesh levels. 

Jameson^ was the first to apply Brandt’s FAS Multigrid scheme to the Euler 
equations, and it is this scheme upon which much of the present work is based. The 
success of Multigrid acceleration for the Euler equations is based on the premise that 
wave propagation, in addition to error damping, determines convergence rates. As the 
problem is solved on increasingly coarser mesh levels, larger time-steps may be used, 
thus propagating the errors out of the domain at a faster rate, and speeding convergence. 
Used in conjunction with Jameson’s hybrid Runge-Kutta schemes, high frequency error 
damping, characteristic of these schemes, is also used. The time-stepping schemes are 
tailored to provide good high frequency damping characteristics, which, when applied 
on the coarse mesh levels, help to eradicate the low frequency errors. 

Multigrid Acceleration has proven to be a robust and reliable tool for the Euler 
equations^ 2,3,6,12-16 ! with recent advances also being made for the Navier-Stokes 
equations.^ 1 3-20 ^ However, improvements are still required, especially for the Navier- 
Stokes equations. The addition of the shear-stress terms and the high mesh stretching 
which accompany these problems tend to decrease the efficiency of Multigrid. The 
theory behind Multigrid for the Euler and Navier-Stokes equations has not reached as 
high a level of development as it has for elliptic problems. It essentially lacks a general, 
theoretical approach, and many times tends to be problem dependent. 

An area which to date has not been explored to any large degree is the effect of the 
coarse-mesh boundary conditions on convergence. Enforcing the boundary conditions 
on the coarse-mesh levels with coarse-mesh accuracy will affect the fine-mesh solution 
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due to the higher truncation error on the coarse levels. For interior points, forcing terms 
are added to the equations for this very reason, and it is the present hypothesis that 
similar terms are required for the boundary conditions. 

It is therefore the objective of this work to obtain a correct, consistent formulation 
of the coarse-mesh boundary conditions and to determine the effect of these boundary 
conditions on the convergence of both inviscid and viscous problems. As an initial 
step in obtaining these conditions, a new, general approach to formulating the Multigrid 
equations is given. The equations on the coarse mesh levels are viewed as a filtered 
sub-set of the fine mesh equations, since certain information is, in effect, filtered out on 
the coarse levels due to the lack of mesh resolution. In this context, a filtering operator 
is first defined. The coarse mesh equations in partial differential form are then derived 
by filtering the original partial differential equations one or more times. The 
specification of a discrete filter then gives the procedure for transferring the fine-mesh 
solution to the coarse mesh. The coarse-mesh equations in discrete form are 
constructed through a finite-volume approximation of the filtered coarse-mesh 
equations. Although this method may be used to obtain any number of coarse-mesh 
discretizations based on how the filter is defined, one choice reduces the present 
formulation to the standard Jameson Multigrid scheme. Further, the present analysis 
leads to the correct formulation of boundary conditions on the coarse levels, without the 
coarse-mesh truncation error affecting the fine-mesh solution. 

The work presented in this thesis is divided into the following chapters. In Chapter 
2, the governing equations are presented. The 2-dimensional Euler, and full, laminar, 


6 


Navier-Stokes equations are described along with their non-dimensionalization with 
respect to ffeestream conditions. The equations are then given in a general, non- 
orthogonal coordinate system. Chapter 3 gives the important aspects of the numerical 
method. Details of the finite-volume cell-centered spatial discretization, 4-stage hybrid 
Runge-Kutta time-stepping scheme, and blended 2nd and 4th difference artificial 
dissipation model are described. In Chapter 4, the boundary conditions for inviscid and 
viscous flow are given. Both a characteristic variable formulation and a Riemann 
Invariant formulation are given for the far-field boundaries. Chapter 5 presents the 
methods that are used to accelerate convergence. Jameson’s Multigrid scheme is 
presented, and implicit residual smoothing for 2-dimensional problems is described. 
The general Multigrid formulation followed by the derivation of the boundary transfer 
operator and the coarse-mesh boundary conditions is given in Chapter 6. In Chapter 7, 
results are presented and discussed for inviscid channel flow over a bump for subsonic, 
transonic, and supersonic speeds, flow over a viscous, flat plate, viscous subsonic 
channel flow over a bump, and inviscid and viscous flow over a VKI gas turbine rotor 
blade. Chapter 8 summarizes the work that was performed and presents the 


conclusions that were drawn. 
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CHAPTER 2 - GOVERNING EQUATIONS 


2.1 Navier Stokes Equations 

The two-dimensional, unsteady, compressible Navier Stokes equations may be 
written in conservative form as follows: 


au 3F 3G 3R 3S 

dt dx dy dx dy 


(2.1) 


where : 
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P u 


pv 

u = < 

pu 

► F = • 

pu z +p 

* G = 1 

puv 

pv 

1 

puv 


pv z +p 


E 

w J 


(E + p)u 

v 4 
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UT XX + Vt xy - q, 

t -J 


UT Xy +VTyy -q y 


( 2 . 2 ) 


and where, based on Stokes Hypothesis: 


2 _ du dv 

Cxx ~ 3 ^ 2 dx dy ) 

. du dv . 

^ q " = ' K 


2 dv du . 

x yy ~ 3 ^ (2 d y dx } 

dT 


(2.3) 


dx 


q y = - K 


§T 

dy 


with density, p , cartesian velocity components, u , v , total energy per unit volume, 
E , pressure, p , temperature, T , viscosity coefficient, |i , and the coefficient of 
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thermal conductivity, k . Pressure is defined by the equation of state for an ideal gas : 


p = pRT = (y - 1) 



yp(u 2 +v 2 ) 


(2.4) 


where y is the ratio of specific heats. Assuming laminar flow, the equations are closed 
by using Sutherland’s Viscosity Law relating the viscosity to the local temperature^ 21 ' 


P = 


3 



T + C2 


where for air and moderate temperatures: 


(2.5) 


c 1 = 1.458X10 -6 kg/(m sec °K '*) C 2 = 1 10.4 °K 

The coefficient of thermal conductivity is expressed in terms of the local viscosity and 
the Prandtl Number P r : 


K = 


CpP 

Pr 


( 2 . 6 ) 


where Cp is the specific heat at constant pressure. 

In this work, the governing equations are non-dimensionalized with respect to a 
reference length, L, and freestream conditions, p« , V„ , T„ , and . 


x 





t = 


tV„ 

L 


* 

P 


P 

PooV 2 





= _p_ 
p~ 


p 
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Where * denotes non-dimensional quantities. Under this non-dimensionalization, the 
Euler equations retain their original form and a constant, 1/ ReL, appears before the 
viscous terms in the Navier-Stokes equations. 

au* a^_ 9G^_ _ _1_ 

at’ + ax* 3y* Rc l 

where ReL is the reference Reynolds number defined as: 


8R* as* 

dx* dy* 


(2.7) 


Re L = 


P- u„ L 
P~ 


(2.8) 


The non-dimensional temperature is derived from the equation of state (Eq. 2.4) : 


7 Mi p* 

* * 

P 

And the non-dimensional viscosity and coefficient of thermal 
derived as: 


(2.9) 

conductivity may be 



1 + C 2 / T,. 
T + C 2 /T 00 


K 


P 

P r (y- 1) M i 


( 2 . 10 ) 


where T«, is a free-stream reference temperature. The equations for the non- 
dimensional pressure and other thermodynamic variables remain the same as their 
dimensional counterparts. For convenience, the * is dropped from the notation at this 
point, and all variables can be assumed to be non-dimensional. 


In order to perform computations over general mesh topologies, the above cartesian 
(x,y) system of equations can be transformed to a non-orthogonal system ( \ , q ), 
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where: 


5 = £(x,y) n = Tl(x,y) 


Equation (2.1) becomes: 


<• ”1 

u 

+ 

d 

' ?,F + 5,G ' 

+ ± 

ti x F+n y G 

J 

„ J 



J 

b 4 

*n 

j 

b 4 


1 

a 

' 4xR + ^yS ' 

+ d 

q x R + rj y S 


Re L 


J 

w d 

an 

J 

b J 



where J is the Jacobian of the transformation and is given by: 

j = h yn- x nn 

Given the metric relationships: 


\ *>y ' 


x $ y$ 

*1x Tly 


n * 


the Navier Stokes equations in general coordinates become: 


Where: 


dU_ + 8F + dG^ = _1_ f dR_ as/ 

at as an Rc l [ a^ + an 


U = U/J 

F = (^F + ^ y G)/J = y T1 F-x 11 G 
G = ( T| X F + q y G )/ J = x^G - y^F 
R = ( £ X R + )/ J = yt|R — 

S = ( r| x R + t| y S )/ J = x^S - y^R 


( 2 . 11 ) 


( 2 . 12 ) 


(2.13) 


(2.14) 


(2.15) 
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CHAPTER 3 - NUMERICAL METHOD 
3.1 Finite Volume Formulation 

A discrete form of Eq. (2. 14) is found using a control volume formulation in which 
the solution values are stored at the cell centers. First, Eq. (2.14) is put into integral 
form for the computational control volume, V: 

Ajj. UdV + JJ. (F*+G„)dV»JJ.(R5 + Su)dV (3.1) 

Then, integration over the volume, A^Arj , gives : 

^ jj ■ U dV + J [ F(5 + A^) - F($) ] dri + { [ G(il + At)) - G(r|) ] d^ = (3.2) 

= J [ R(5 + A^) - R(£) ] dr| + j [ S(T] + At)) - S(T1) 1 

A A 

Applying the Mean Value Theorem to the control volume in Figure 3.1, F and G for 
each side of the control volume may be given by their values at each particular cell face 
center. Eq. (3.2) becomes: 

JLjJ. UdV+[F i + , *.j -Fj.vi.j] ATI + [G i>j + V 4 -Gi.j-viJA^a 

= [ Rj + V4 , j — Ri - V4 . j 1 AT| + [ S^j + 1,4 — . j - 14 ] A^ 


(3.3) 
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2 

i + , j 

4 iJ-'A 1 

Figure 3.1. Cell Centered Control Volume Notation 

Integrating the first term in Eq. (3.3) produces an ordinary derivative in time for the 
solution U at each cell, (i j): 

A|J.UdV=|.J/U<r'd5dT,)=[AV^j (3.4) 

l Ji. j 

where AV is the cell area, which is defined in this work by the cross-product rule as : 

AV = -1/2 [ (x 2 - X 4 >(yi - y 3 ) - (xj - x 3 )(y 2 - y 4 ) ] (3.5) 

The Jacobian of the general coordinate transformation describes the ratio of the cell 
areas between the ( % , tj ) coordinate system and the ( x , y ) coordinate system. 
Since the cell area in the transformed computational domain is equal to 1 

( A£ = Ar| = 1 ) , AV can also be defined as the inverse of the Jacobian : 

AV = y = y T! - yj, Xt, (3.6) 

The second term in Eq. (3.3) may be expressed in F and G as: 

+ At) = [( y^F - x^G )j + ^ , j - ( y^F - x^G X _ ^ , j] At| = (3.7) 

= ( FAy - GAx )j + v* . j - ( FAy - GAx ) 4 _ i* , j 

Treating the remaining terms of Eq. (3.3) in the same manner produces a spatially 
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discrete system of ordinary differential equations for the cell shown in Figure 3.1 : 


dU 

dt J i>j 


— i-[(FAy-GAx) i + *, j -(FAy-GAx) i _ 1/iij + (3.8) 

A Vi t j 

+ ( GAx - FAy )j , j + 1/4 - ( GAx - FAy ) ( , j _ * + 
-(RAy-SAx )i + v 4 ,j + ( RAy-SAx )| _ va . j + 

- ( SAx - RAy )i , j + + ( SAx - RAy )i , j _ ^ ] 


where: 


Ax i + l/4 ,j =X 2 -X! 
AXj , j + i A = X 2 - X 3 
Axj _ Vi . J = X 3 - x 4 
Ax; t j _ — Xj — x 4 


Ay i+ V4.j = y2 -yi 
Ayi. J + ^ = y2-y3 
Ayj - V4 . j = y3 -y4 
Ayi.j-v4 = yi -y4 


(3.9) 


In the present work, the convective terms, F* + y* t j and Gj + y* < j , are obtained from the 
average of the conservative variables, U, existing in the two cells adjacent to the face. 
For linear problems, this is equivalent to averaging F and G. However, since F and G 
are non-linear functions of U, Turkel^ claims that averaging U instead helps to couple 
the even and odd points throughout the domain. To save on the number of 
computations that must be performed at each point, pressure is stored at the cell centers 
and averaged to the face. 


(pu)i + V4,j 


(pur 


+ Pi + V 4 , j 


Fi + V 4 , j = 


i + l /2 , j 


(pu)(pv) 
p 

Ei + 1/4 , j 


i + Vi , j 


r Gi + V4 . j = 


(Pv)i + V 4 .j 


(pu)(pv) 

P 


(pv) 2 


i + Vi , j 
+ Pi + Vi , j 


(3.10) 


i + Vi . j 
Ei + vi , j 
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where 

(p)i + vi . j = y [(P)i+i,j +(P)i.j ] (3.11) 

(pu) i + ^.j = Y [(pu)i+i.j +(pu)i,j J 

(pv) i + Vi . j = y [ (P v )i + 1 , j + (P v )i . j j 
(E) i + *,j = y [(E) i + 1<J + (E) ii j J 

(p)i + >/4 . j = y (P)i + 1 . j (P)i . j j 

Quantities at the other faces are defined in the similar manner. 

The evaluation of the viscous fluxes for each cell requires that the first differences 
of u , v , and T be defined at the center of each cell face. Spatial discretization of 
these derivatives for the point ( i + x h , j ) is performed by integrating over a surface 
bounded by the two adjacent cell centers, (i,j) and (i+1 j), and the endpoints of their 
dividing face, [17) shown in Figure 3.2. 

3 

7 |\ 

/ \ 

4 < » > 2 

\ / 

I 

i 

Figure 3.2. Control Volume for Evaluation of Shear Stress Derivatives at face (i+l/2j) 

For this control volume formulation, the following relationship may be derived for any 
variable, f : 
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JJ.fxdV.J^fdy (3.12) 

Jiv f y dv= -; w fdx 

Where V and 9V are the surface area and boundary, respectively, of the control 
volume in Figure 3.2. The line integrals are approximated by trapezoidal integration, 
which gives a discrete representation for the first derivatives of f: 


df 

9x 

9f 

By 


1 Z (fk+i +f k )(yk + i -yv ) 


2 AV 


(3.13) 


k=l 

1 4 


2 AV 


Z ( fk+ 1 + fk ) ( x k+ 1 ~ x k ) 


k=l 


Where 


fj = ft x 5 = x i Y5 = Yi 

Using this method, the shear stress terms can be evaluated in their cartesian form 
without the need for a general coordinate transformation. When evaluating these 
derivatives at boundaries, such as the (i,j-l/2) face of the first interior cells, the center of 
the image cell must have specified coordinates for the above finite-volume formulation. 
To avoid difficulties of defining these coordinates, such as in areas of high curvature, a 
different approach is utilized at the boundaries. The viscous control volume of the 
(ij-1/2) face of a cell along the boundary is collapsed to a triangular volume, shown in 
Figure 3.3, in which the first derivatives required at the boundary are obtained directly 
by integration over the truncated control volume. This results in a lst-order accurate 
evaluation of the boundary shear stress terms, which is essentially equivalent to 
performing one-sided differencing there. 
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Figure 3.3. Control Volume For Boundary Shear Stress Evaluation 


The finite-volume, cell centered formulation for the Navier-Stokes equations 
defined in this section is equivalent to a central difference scheme, and is 2nd order 
accurate for sufficiently smooth meshes J 6 ' 

32 Time Stepping Scheme 

The time-marching algorithm used in this research is a 4-stage Runge-Kutta 
scheme. Typically, multi-stage time-stepping schemes such as these are chosen for 
their high order of time accuracy. However, since the goal of this research is steady- 
state computations, time accuracy is not a requirement. 

This 4-stage scheme is one of a class of hybrid multi-stage schemes all of which 
were developed by Jameson 16 ' specifically for their damping and stability 
characteristics. The coefficients of the scheme are tailored to give the maximum 
allowable time step and to provide optimum damping of the high frequency error 
modes. The attenuation of these modes is essential for the success of rapid steady-state 
convergence. 


The semi-discrete system of ordinary differential equations in Eq. (3.3) may be 
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rewritten as 


^ = -1[ C (U) + V (U) - D (U)] (3.14) 

dt AV 

Where C (U) is the convective flux, 

C (U) = (FAy - GAx)j + * . j - (FAy - GAx)i _ * . j + (3.15) 

+ (GAx - FAyX , j + v* - (GAx - FAy)j , j _ * 

V (U) is the viscous flux, 

V (U) = - (RAy - SAx)i + * , } + (RAy - SAx)j _ j - (3-16) 

- (SAx - RAy); t j + % + (SAx - RAy); , j _ v* 


and D (U) is the artificial dissipation which will be defined later in section 3.3. 


The 4-stage Runge Kutta scheme used to complete the discretization of Eq. (3.11) is 
implemented in the following manner : 


u (0) = u n 


u (1) = u< 0) “ lAt 

AV 

C (U (0) ) + V (U (0) ) - D (U (0) )] 

tr(2) _ y(0) a ^ At 

AV 

C (U (1) ) + V (U (0) ) - D (U (0) )] 

U(3) = U (0) tt 3 At 

AV 

C (U (2) ) + V (U (0) ) - D (U (0) ) j 

u< 4, = u«» 

AV 

C (U (3) ) + V (U (0) ) - D (U (0) )] 

u n+ 1 =u (4) 



The superscripts in parentheses refer to the particular stages of the scheme, and where 
otj are the coefficients designed to give optimum stability and damping to this 
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scheme.^ 


a! = 1/4 a 2 = 1/3 (X 3 = 1/2 04 = 1 

Note that the viscous fluxes and the artificial dissipation are evaluated once and held 
constant throughout the time step. Other Runge-Kutta schemes may possess better 
damping characteristics or allow larger time steps, but the single evaluation of the 
artificial dissipation in this 4- stage scheme produces substantial savings in 
computational time. 

This 4-stage Runge-Kutta scheme has been used with good success in many 
steady-state computations, and it’s high frequency damping 
characteristics make it ideal when used in conjunction with muitigrid acceleration.^ 


3.2.1 Stability Criteria 

Turkel^ derives the stability limit for the Euler Equations based on the maximum 
eigenvalues of the convective Jacobian matrices. Ignoring the effect of artificial 
dissipation on the equations, the maximum allowable time step for a general, multi- 
stage scheme is given there as : 


At £ 

luyT, -vx^l + |vx^-uy$| +( 

An approximation of which is : 


AV (CFL) 

A + 4 + A + 4 + 2 1 x^ Xt1 + y § y n I )*a 


(3.18) 


At ^ AV(CFL) 


(3.19) 


Where a is the speed of sound, AV is the cell area, and CFL is the Courant number for 
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stability based on a Fourier analysis of the 1-D model equation : 


Ut ^xxxx — 0 (3.20) 

and Xyy are the maximum local wave-speeds in the % and r| directions, 
respectively, and are defined as : 




[l uyr, - VXr, | 
[l VX^ -uy^l 


+ a Vyn + x n ] 
+ aVy| + 4) 


(3.21) 


The maximum CFL number that may be used is dependent on the type of hybrid 
Runge-Kutta scheme that is used. Factors affecting this stability condition are the 
Runge-Kutta coefficients, cq, the number of evaluations of the artificial dissipation, the 
number of Runge-Kutta stages, and the amount of smoothing, p . The maximum CFL 
number for an m-stage scheme in which the coefficients are optimized for the largest 
time step is^ : 


CFL < m-1 (3.22) 

For the present 4-stage scheme with a single evaluation of the artificial dissipation, Eq. 
(3.17), with p. = 1/32, gives the condition : 

CFL < 2.6 (3.23) 


For the Navier Stokes equations, diffusive effects must be considered in the stability 
analysis as well as the convective effects. In Reference [18], the Navier-Stokes time 


step is given as : 
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At = CFL 


Ate At v 
At<; + At v 


(3.24) 


where Ate is the convective time step given in Eq. (3.15), and At v is the viscous time 
step : 


At v =K v y- (3.25) 

K v is a specified empirical constant that weights the importance of the viscous terms 
over the convective terms and is given as 0.25. [18] Xy is the sum of the maximum 
spectral radii in the % and T| directions of the viscous operator in the Navier-Stokes 
equations : 


Where 


Xy — Xy^ + Xy^ 


(3.26) 


\ - 
= 


1 


Re L AV 

1 

Re L AV 


-Hi-Al 2 + f-AlAm 
P r p 3p 

Am 2 + -^-AlAm 
P r p 3p 


(3.27) 


and where A1 and Am are the lengths of the cell in the £ and r\ directions, 
respectively: 


ai = Vh + n 
Am = v*; + yj, 


(3.28) 
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For the viscous cases considered in the present work, the maximum CFL number 
used with the above viscous time step criteria was approximately the same as those for 
the inviscid cases. 

3J Artificial Dissipation 

Central difference schemes require the addition of artificial dissipation to damp the 
high frequency error modes that occur as a result of the odd-even point decoupling. 
Artificial dissipation is also required for shock capturing and to damp oscillations in 
other high pressure gradient regions, such as stagnation points. Even though the 
viscous terms in the Navier-Stokes equations provide physical dissipation and would be 
capable of resolving the structure of a shock, mesh spacing in the region would need to 
be on the order of the molecular mean-free-path. This is not practical from a 
computational standpoint, and therefore artificial dissipation is still required to capture 
shocks. Also, artificial dissipation is required to damp instabilities that may occur in 
regions dominated by the convective terms. 

In the present work, the dissipation model is based on that introduced by Jameson, 
Schmidt, and Turkel. [5] It employs modifications by Swanson and Turkel, [221 and 
Martinelli, t23] to improve accuracy and to increase convergence rates for viscous 
solutions. It is a blended 2nd and 4th difference adaptive dissipation scheme that 
provides a base level of 3rd order dissipation throughout smooth regions of the flow, 
and decreases to 1st order in the vicinity of shocks and other high pressure gradient 


regions. 
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The conservative form of the equations is maintained by evaluating the dissipation 
at each cell center by summing the dissipative fluxes in each coordinate direction. The 
basic form of the dissipation, Dj , j , for each cell is: 


D i.j= d 5 _d l +d i| ~ d Ij] u i.j 


(3.29) 


Where Uj , j is the vector containing the conservative variables, and d 2 and d 4 are 
the 2nd and 4th adaptive difference operators. For cell ij, shown in Figure 3.1, these 
are defined as : 


d^U = V^ 

<*;u=vJ 


^i + Vi,j £ i + 'h , j J j f j J 

^ i + >/i . j e i + >/4 , j ] V^U i , j J 

^ i * j + V* + > 

j 

^ i , j + x h £ i , j + x h 


(3.30) 


Where A and V are the standard forward and backward finite differences, 
respectively. 


A^U = U i + 1 . J -U i . J V^U = U i>J -U i _ 1>j (3.31) 

A T1 U = U i . J + 1 -U i , J V n U = U i . j -U i . J _ 1 


Xj + vi , j and Xj , j + v* are local variable scaling factors averaged to the cell face 
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i + V4 . j 
£ i . j + Vi 




+ S..i] 


(3.32) 


In the original model, 



AV 

At* 


where At* is the inviscid time step due to a unit CFL number, or. 


(3.33) 





(3.34) 


where X.^ and X^ are the maximum local wave-speeds in the t, and T| directions 
defined in Eq. (3.21). The adaptive coefficients, e (2) and e (4) control the blending 
of the 2nd and 4th differences in the dissipative operator : 


<» -r<» 

e i + vs . j = K max 


Vi-i 


* • j 


V 


i + 1 , 



(4) 

£ i + Vi , j 


= max < 


0 , 



- e 


( 2 ) 

i + Vi , j 


> 


J 


e*j + Vi = K<I> max (vi.j-t , Vi,j , v itj+1 
e T! j + V4 = max jo, [ K (4> - e f’ j + j 


(3.35) 


Where v is the norm of the centered 2nd difference of pressure used to locate large 
pressure gradients and turn on the second difference dissipation. 
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< v i.j>5 = 


P i + 1 , j 2P j t j + P i _ j _ j 
P i + 1 , j "*■ 2P i , j + P i - l , j 

P i.j+1 ~ 2P i,j + P i.j-1 
P i.j+l+ 2P i,j+ P i.j-l 


(3.36) 


When the second differences are strong, the fourth differences are turned off by e * to 
sharply resolve shocks. K™ and K <4> are user specified constants which typically hold 


values of: 


K 

K 


( 2 ) 


(*) 


1 

4 1 
1 

256 


2 


to 


50 


The values that are used are based on the best convergence rate that may be obtained 
(highest possible K ’ ) while still maintaining accurate solutions, and on the ability to 
cleanly capture shocks without too much smearing (K™). 


The difference operators in Eq. (3.19) arc applied in two steps. A sweep is made 
through the domain in the £ direction, taking centered first differences of U at each 
cell face for d 2 , and 3rd differences at each face for d 4 . Then another sweep is made, 
taking centered differences of the 1st and 3rd differences, yielding the desired 2nd and 
4th differences, respectively, at cell centers. The same procedure is repeated in the rj 
direction. 


3J.1 Modifications 

Based on an analysis of Eq. (3.25), Swanson and Turkel* 22 * determined that the 4th 
difference operator produces dispersion as well as the required dissipation. If the 4th 
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differences are applied as a sequence of two 2nd differences instead of a 1st and a 3rd 
dif ference, the operator d 4 will produce only dissipative terms. 


D^U = 


(*i.i eHj] 

O 

j * 

C 

II 

’^n]> 

(li.i 


Note that e^ 4) and X are now evaluated at the cell centers, not at the faces. 

The second area of modification is in the evaluation of the variable scaling factor, 

X . In viscous calculations, severe mesh stretching can result from the need to resolve 
the boundary layer. Cell aspect ratios may vary from Al/Am»l near a solid; 
boundary, to 1 in the far-field. This large variation slows convergence, and diminishes 
the accuracy of steady state solutions. Also, these problems increase in magnitude for 
multigrid calculations due to the large difference of the high-frequency modes in the 
two coordinate directions.* 22 ' 

To overcome these difficulties, a new X based on cell aspect ratio was introduced 
by Martinelli* 23 ' It is essentially a combination of the scaling factor in the original 
model with an anisotropic scaling factor suggested by some researchers.* 13 ** 24 ' In the 
£ direction, 


(>n,j \ = h 



(3.38) 


And in the r\ direction, 



Where 


0 < a i£ 1 (3.40) 

If a is 1, then X reduces to the scaling in the original dissipation model given in Eq. 
(3.30). If a is 0, X reduces to scaling in one direction only at each given face. This 
scaling is not recommended as it may decrease convergence rates and create problems 
for multigrid applications 1221 . A value of a that has produced good results 12211181 is 

a = 2 / 3 

For the cases in the present work, solution convergence did not seem to change 
significantly with changing a . 

In the present cell-centered scheme, image cells are included around the physical 
domain to allow the same algorithm to be used for all interior cells. In order to 
determine the 4th difference dissipation terms in cells adjacent to boundaries as in 
Figure 3.4, 2nd differences must be given at the image cells. 

(VAVAU) 2 = (VAU) 3 - 2(VAU) 2 + (VAU)! (3.41) 

Where V and A are the difference operators given in Eq. (3.27), and where : 


(VAU)i = U 2 - 2U, +U 0 


(3.42) 
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0 


Figure 3.4. Cell Notation for Boundary Dissipative Flux 



Since these differences can not actually be calculated ( Uo is not defined) they must be 
assigned in a manner that maintains conservation form and doesn’t introduce instability 
at the boundaries. 

Swanson and Turkel^ 22 ' provide a detailed analysis of various treatments of these 
terms near boundaries. For far-field boundaries, the standard approach is to set the 2nd 
differences in the image cells to zero. 

(VAU)j=0 (3.43) 

A typical treatment at solid boundaries is to set the surface dissipative flux to zero 

(VAU)i = (VAU) 2 (3-44) 

and then either : 

1) (VAU) 2 =0 (3.45) 

or 

2) (VAVAU) 3 = 0 
(VAU) 3 =(VAU) 2 =0 

Option 2 results in zero dissipation normal to the boundary for the first two interior 


cells, 2 and 3. 
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In the present work, both methods were used with equal success, and without 
noticeable differences in the results. 
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CHAPTER 4 * BOUNDARY CONDITIONS 

4.1 Inflow/Outflow Boundaries 

For the channel flow cases, where the outflow conditions are essentially freestream, 
a characteristic formulation of the boundary conditions is used. For the turbine cascade 
case, in order to explicitly satisfy constant total temperature, total pressure, and flow 
angle, a Riemann Invariant formulation is used at the upstream boundary, while a 
characteristic formulation is retained at the outflow boundary. 

4.1.1 Characteristic Formulation 

This approach utilizes a transformed system of equations in which the boundary 
cells are updated by specified and extrapolated characteristic variables. Rewriting the 
Euler equations in a local coordinate frame that is normal and tangential to the 
boundary gives : 

*J. + iF + iO=0 (4.1) 

9t 9n 9s 

Where n and s denote the normal and tangential directions. At these boundaries, it is 
assumed that variations of the flow conditions in the tangential direction are negligible. 


and so Eq. (4.1) can be written as : 
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au dF 
at an 


(4.2) 


The above one-dimensional system of coupled differential equations can be linearized 
and transformed to an uncoupled linear system of partial differential equations: 


aw . - aw 
at an 


=o 


(4.3) 


where W is the vector of characteristic variables, and X is the diagonal matrix of 
eigenvalues. Denoting and q s as the normal and tangential velocities at the 
boundary, where the normal is pointing into the solution domain, X and W are defined 
as : 


q» 0 o 0 1 P P/> 

x = o q» o o w= 4 __ 

0 0 q„+a 0 (1/ 'V2)(q n + p/pa ) 

0 0 o Qn~a r 

L J |0/V2)(-q n + p/pa ) 

The barred terms represent variables evaluated at the linearized state, which in this 
work are taken as the values from the previous Runge-Kutta stage. 

Along the characteristics defined by slope dt/dn= 1/ A^, W* is constant and Eq. 
(4.3) reduces to an ordinary differential equation. If the slope of the characteristic, 
1/Xj, is such that it originates from outside the computational domain, Wj must be 
specified. If the characteristic originates from within the computational domain, the 
value of Wj at the boundary for the new time level, t + At is found by tracing the 
characteristic back to its intersection point with the previous time level, t. Evaluating 
W, at this intersection point will give the new value of W ; since W; is constant along 
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these characteristic lines defined by dt/dn=l/A,. However, in steady-state 
computations, it is common practice to approximate the value of W ; which lies along 
the characteristic by the value contained in the first interior cell (Figure 4.1). 



Figure 4.1. Approximation To W ; For The Characteristics Originating Within The 
Domain 


The boundary conditions are enforced at the boundary, or cell face, instead of the image 
cell (Figure 4.2) to accommodate a consistent formulation of the boundary forcing 
terms that will be given later in Chapter 6 . The conservation variable in the image cells 
are then updated by linear extrapolation from the interior as: 

Ui=2U b -U 2 (4.5) 



Figure 4.2. Boundary Cell Notation 


For subsonic inflow, and an inward pointing normal, n, , \i , and X 3 are 
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positive, meaning the characteristics originate from outside the domain (see Figure 4.3). 
Therefore, W t , W 2 , and W 3 must be specified, and W 4 must be extrapolated from 
the interior. The following equations result: 


Pb P~ 

Pb - -j: - P~ - ~jr 
a a 


(4.6) 


% =£ ls- 

Pb P~ 

+ — =qn_ +3^ 

P a pa 

, Pb Pex 

P a pa 


where the subscripts b, <», and ex refer to the boundary, freestream, and extrapolated 
values, respectively. Solving Eq. (4.6) for the primative variables at the boundary 
gives: 


q S(1 = q s _ (4.7) 

Pb = y [ Pex + P~ + P a ( q„_ - q,^ )] 

Pb = p- + ( Pb - P~ V a 2 
= qn_ + ( P~ - Pb V P a 


For subsonic outflow, with an inward pointing normal, n, only X 3 is positive 
(Figure 4.3), and so W 3 is the only characteristic variable that must be specified. 
Following the same procedure as above, the updated values are given by : 
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q Sb = q Su ( 4 -8) 

Pb = yt Pout + Pex + P a ( qiw ~ qn,, )1 

Pb = Pex + ( Pb ' Pex V 3 
qn b = q^ + ( Pb " Pex V p a 

where the subscript out refers to the exit conditions. 

For the inviscid channel flow cases, p out and are assumed to be at ffeestream 
conditions. For viscous problems, where the boundary layer (or viscous wake) 
intersects the outflow boundary, the characteristic relations are incorrect due to the 
non-isentropic conditions existing there. For these cases, the exit pressure is specified 
and the other variable are extrapolated. 

For supersonic inflow, all four Xj are positive (Figure 4.3), and therefore all four 
W, must be specified, which reduces to specifying the primative variables: 

q Sb = q s _ ( 4 - 9 ) 

Pb = P~ 

Pb = P~ 
qn b = qn_ 

For supersonic outflow (Figure 4.3), all X* are negative, and so all Wj must be 
extrapolated: 

= qex 
Pb = Pex 
Pb = Pex 

qk = 


(4.10) 
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t 



Figure 4.3. Characteristics 


t 




At Far-Field Boundaries 


4.1.2 Riemann Invariant Formulation 

In this formulation, an upstream-running Riemann invariant is used in conjunction 
with specified total pressure, total temperature, and flow angle to solve for the updated 
inlet conditions.' 25 ' The Riemann invariant, R - , is based on total velocity, Q, and is 
approximated by extrapolation from the first interior cell : 


R“=Qex 




Where a is the local speed of sound 


a = 



(4.11) 


(4.12) 


Computing R in this fashion represents extrapolation along streamlines versus normal 
to the boundary. Based on isentropic relations and the specified total temperature at the 
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inlet, T 0 _ , the inlet speed of sound can also be given by : 


a b = 


( Y~ 1 )C p T 0 _ - Qi( *^r— ) 


V* 


(4.13) 


Where Cp is the non-dimensional specific heat 

1 


Cp = 


Mi(y-l) 


(4.14) 


and Qb is the unknown inlet total velocity. Using Eq. (4.12) and Eq. (4.14), a 
quadratic in Qb results. Being an inflow boundary, the negative root is non-physical, 
and therefore Qb is given by the positive root: 


Qb = 


(Y-1)R- + 


- \2 


4( y + 1 )CpT 0 _ — 2( y— 1 )( R _ ) 


y+1 


(4.15) 


Pressure and density at the inlet are then calculated using isentropic relations and the 
specified inlet total pressure : 

l-Y/< Y- 1 ) 


Pb — Po_ 


Pb — Po. 


1 + 


Pb 

Po. 


(Y- 1 ) 
2 

l/Y 


r •N 

Qb^ 

2 

a b 

k. J 



(4.16) 


The cartesian velocity components are then given by Qb and the specified inlet flow 
angle, pi : 
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u = Q b cosp 1 (4.17) 

v = Q b sinp 1 


At the outflow, the characteristic formulation (section 4.1.1) is used with a specified 
exit pressure. The exit pressure in the turbine case is h ased on a given exit isentropic 
mach number M out]i , assuming that p 0ou = p 0< _ : 


Pout — Po. 




-Y/(y-n 


(4.18) 


4.2 Solid Wall Boundaries 

For a solid wall in viscous flow, the physical boundary conditions required to close 
the system of governing equations are zero velocity and either specified surface 
temperature or heat flux: 


Q = 0 

T _ T 3T 

T — i wall or — — 


(4.19) 


9n 


9T 

an 


I wail 


where Q is the total velocity and n is the normal to the surface. The physical boundary 
conditions required for inviscid flow is flow tangency at the surface (zero flux through 
the surface): 


Q-n = 0 (4.20) 

For a numerical scheme, however, other boundary conditions are required to close the 
discrete equations, and these must be specified in a manner consistent with the physics 
that are occurring at the boundaries. The implementation of these numerical boundary 
conditions at a solid wall is given below for inviscid and viscous surfaces. 
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For a cell-centered, finite volume scheme, the only contribution to the momentum 
equations at a solid wall is the pressure. The pressure may be obtained by extrapolation 
from the interior based on the pressure derivative normal to the boundary 151 : 


ap ( x^ + y^ )|| + p( y^u - x^v )( x^v - y^u ) 

(x^ + y^) 

giving 


(4.21) 


Pb = P2 


ap 

an 


(4.22) 


Where the subscript b denotes the boundary point and the subscript 2 denotes the first 
interior cell as before (Figure 4.2). For inviscid flow, flow tangency is specified by 
setting 


q„.=0 (4.23) 

qs b = Qs* 

where q„ and q s are the normal and tangential velocity components. Density is 
obtained by isentropic relations based on the assumption that the entropy gradient 
normal to the boundary is zero: 


Pb = P2 



(4.24) 


For the viscous cases, the normal pressure derivative at the boundary is assumed to be 
zero in the present work. A no-slip condition is specified for the velocity : 
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<k =° 
q** =o 

and a specified wall temperature, Tg, is used to define the density : 

yMipb 

Pb = 


For the flat plate problem, an adiabatic surface is specified: 


9T 

dn 


= 0 -4 T b =T 2 

j b 


(4.25) 


(4.26) 


(4.27) 


For the viscous bump and turbine problems, the temperature at the boundary is set equal 
to the ffeestream total temperature: 


(4.28) 
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CHAPTER 5 - CONVERGENCE ACCELERATION 
5.1 Multigrid Method 

Multigrid was developed by Brandt 110,11,1 to increase convergence rates for the 
solution of steady, elliptic-type problems, and was later applied to the Euler equations 
by Jameson. 1 61 Jameson used Multigrid directly with his hybrid Runge-Kutta schemes 
on all mesh levels, solving the same fine-mesh equations with the addition of forcing 
terms to maintain fine-mesh accuracy. In this chapter, a brief discussion of Brandt’s 
FAS scheme for steady problems will be given. Then his notation will be used to 
describe Jameson’s procedure for applying Multigrid to the time-dependent Euler 
equations. 

A general, non-linear steady equation is given in Brandt’s notation as : 

L h U h = F h (5.1) 

where h denotes relative mesh spacing, L is the discrete spatial operator, and F is a 
forcing term which is usually zero on the finest mesh level. The corresponding 
coarse -grid equation is : 


L2hU2h = P2h 

where the forcing term is based on the transferred fine-mesh residual (Fh - L^Uh): 
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F2h = In’ ( Fh - 4U h ) + l<2h ( Ih’Uh ) (5.3) 

I h is the fine-to-coarse transfer operator (interpolation scheme). Changes in the 
solution are transferred back to the fine mesh to update the solution : 

Uiw - Uh*. + AU a (5.4) 

where 


AU 2h = u 2h - I^Uh^ (5.5) 

and where I 2 h is the coarse-to-fine transfer operator. This transfer introduces high 

frequency errors into the solution on the fine grid, and therefore, it is important that the 

relaxation scheme (or time-stepping algorithm for unsteady equations) possess good 

high- frequency damping characteristics 

For the unsteady Navier-Stokes equations using Runge-Kutta time-stepping, 
Jameson’s Multigrid equations are of the form: 

U 2 h - - 0tAt( I^hUa - Fa ) (5.6) 

where the coarse -mesh forcing term is : 


Fa = ( F h - L h ulT ) + L 2h ( lg h Uh ) 


with F h = 0. The quantities in parentheses refer to the stage of the Runge-Kutta 
scheme, and L represents the discrete spatial operator given in Eq. (3.1 1) : 


L 2 hU 2 h = 

ih ih AV 


C(U 2h ) + V(U 2h )-D(Ua)] 


(5.7) 


The forcing term, formulated in the above fashion, is necessary in order that the fine- 
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mesh solution will not be affected by the coarse-mesh truncation error. If the fine mesh 
converges (F h -L h U h =0), the coarse mesh converges also. Looking at the 1st 
Runge-Kutta stage on the coarse mesh: 

Ujh = u! - aAt( L 2h u‘J - Fa ) (5-8) 

= U» - otAt[ L 2h U; - I? 1 ( F h - L h Uh ) - L 2h ( ifu" ) ] 

= u5-«At[ - I? ( F h - L h U ( h )] 

— i r (0> 

Or, in other words, 3U/ dt = 0 on the coarse mesh. 

For Jameson’s cell-centered scheme, the transfer operator I h is a volume 
weighted average over the 4 fine-mesh cells that make up a given coarse-mesh cell 
(Figure 5.1). 
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Figure 5. 1 . Fine-to-coarse Transfer 


The transfer of the solution and the residual from the fine mesh are therefore given as: 
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Solution transfer: 


I U hi AVh, 

I^ h U h = ^ 

ZAV h , 


Residual transfer: 


(5.9) 


I?( F h - L„U h ] = 


IAV, 

1=1 


- 1 
AV; 


[C(U h ) + V(U h )-D(U h )] 


IAVj 


I Fluxi 

is: 

1=1 


= £ Residualj 

i=i 


(5.10) 


This residual transfer approximates the integral equations on the coarse mesh cell 
(2 A^Atj) with fine mesh accuracy, and maintains their conservative property. 

Bi-linear interpolation is used to transfer the coarse mesh solution changes back to 
the fine mesh (Figure 5.2) : 


All* . j = AU, + -L [ AU 2 + AU 4 ] + AU 3 


(5.11) 


AU i+ , . j = ^ AU 2 + -L [ AU, + AU 3 1 + -1. AU 4 
AU, , j + , = AU 4 + 1 AU i + AU 3 1 + TZ AU 2 


AU 1+1 . J+1 =4-AU 3 + ^[AU 4 + AU 2 ] + -^AU, 
lu lo 16 



Where the subscripts i ,j denote the indices for the fine-mesh cells and 1-4 denote the 
coarse-mesh cells, and where AU is the coarse-mesh solution change given in Eq. (5.5). 
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Figure 5.2. Coarse-to-fine Transfer 


The cycling process employed in this work is referred to as a simple Saw-Tooth 
cycle (Figure 5.3). One Runge-Kutta integration is performed on each mesh level until 
the coarsest mesh is reached, and then the solution changes are transferred back to the 
finest mesh without any integration steps being performed on the way down. This may 
not be the optimum cycling strategy, but it has been found to be computationally 
efficient for a number of flow calculations.' 6 " 13 ' ' 7 ' ' 26 ' 
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Grid Level 
1 



Runge-Kutta integration 

* Fine-to-coarse transfer (I^ 1 ) 
Coarse-to-fine transfer (I^,) 


Figure 5.3. Saw-Tooth Multigrid Procedure 


52 Implicit Residual Smoothing 

The concept of implicit residual smoothing for Runge-Kutta schemes was 
introduced by Jameson^ to permit the use of larger time steps and therefore increase 
steady-state convergence rates. The residuals at each point, R^, are replaced by 
smoothed residuals, Ry by implicitly solving the equation: 

Ri . j = R i . j + (V^)Rj , j + e n , j (5. 12) 

where V and A are the first forward and first backward differences, respectively, that 
were defined in Chapter 3, and e, and are smoothing coefficients in each coordinate 
direction. This could be thought of as solving the partial differential equation: 
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3R _ e t 

a 2 R 

£ 

+ _JL 

3 2 r 

5t At 

J 

At 



(5.13) 


which is similar in form to the 2-D Heat Equation, and possesses the same properties. 
This equation (Eq. 5.13) therefore has the effect of smoothing the residual at each point 
based on all the residuals in the entire domain. The global influence on each residual is 
determined by e, a "diffusivity” term. Discretized implicitly, Eq. (5.13) is the same as 
Eq. (5.12). Rearranging Eq. (5.12) gives: 


(1— )Ri. j - R i, j (5.14) 

of which a factored approximation is : 


( 1 -e ft )( 1 - e n V n A,, )Rj . j = Ri . j (5.15) 

In the present work, Eq. (5.15) is applied as a combination of one-dimensional 
smoothing steps : 


( 1 “^V^A^ )R*. j = R i. j (5.16) 

( 1 — £,, ^T) &r\ )^i . j = . j 

Where R* are the intermediate smoothed residuals after the sweep in the \ direction. 


Neglecting the effect of artificial dissipation, Jameson {6] showed that stability can 
be maintained for any CFL number so long as £ satisfies the requirement : 


e 



(CFL) 2 

(CFL*) 2 


(5.17) 


where CFL* is the stability limit for the unsmoothed scheme. However, using CFL 
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numbers that are too large may severely decrease high frequency error damping and 
reduce convergence rates. 171 Also, with non-linear equations containing artificial 
dissipation, other factors may enter in and create instabilities. Jameson 161 suggests that 
optimum convergence rates are usually obtained for 

CFL = 3(CFL*) (5.18) 

For even-stage hybrid Runge-Kutta schemes, stability is maintained when the 
smoothing is applied during every even-numbered stage, while for an odd- stage 
scheme, it is applied during odd-numbered stages. 

Locally varying and may be used at each point to optimize the smoothing for 
highly stretched meshes. 11811191 However, in this work, good results were obtained by 
using constant coefficient values of : e = e = 2.0 



CHAPTER 6 - MULTIGRID BOUNDARY CONDITION FORMULATION 


In Multigrid applications, a topic that has received little attention is the effect of the 
coarse-mesh boundary conditions and boundary transfer operators on the accuracy and 
convergence of the solution. In fact, the procedure for dealing with these issues is 
rarely mentioned, and if so, not in a very detailed fashion. Merely enforcing the fine- 
mesh boundary conditions on the coarse mesh levels will introduce coarse-mesh 
truncation error into the fine-mesh solution. Therefore, a consistent formulation is 
required for the boundary conditions just as it is for the interior Multigrid equations. 

In this section a new way of viewing the Multigrid acceleration process will be 
introduced in order to derive a general approach to obtaining the coarse-mesh Multigrid 
equations. This, in tum, will lead to a consistent formulation of the coarse-mesh 
boundary conditions. Rather than formulating the Multigrid acceleration scheme 
directly from the discrete problem, the governing partial differential equations used in 
the Multigrid acceleration scheme will first be constructed. This is accomplished by 
introducing a filter operator which is applied one or more times to the original system 
of partial differential equations. The resulting series of filtered partial differential 
equations may then be discretized leading to a consistent series of coarse mesh 
equations to be solved in the Multigrid acceleration scheme. This approach is 
motivated by the observation that both the coarse mesh solutions and the coarse mesh 
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discrete equations solved in a Multigrid acceleration scheme are in fact filtered 
approximations to the corresponding fine mesh solution and discrete equations, 
although they are rarely viewed in this manner. The advantage of the present approach 
is that the filtering process is clearly defined through the choice of the filter operator. 
Further, since in the present formulation a general filter is initially assumed, a wide 
class of different multigrid discrete formulations are possible by simply changing the 
actual filter that is used. Finally, through a clear definition of the filtering process 
between mesh levels, the proper formulation of the corresponding coarse mesh 
boundary conditions follows directly. 

In the discussion which follows the filtered form of the governing equations will be 
derived through the introduction of a general filter operator. It will then be shown that 
with one possible choice of this filter the corresponding discrete equations reduce to the 
FAS Multi-grid scheme proposed by Jameson. The remainder of this chapter will then 
center on the proper formulation of Multigrid boundary conditions, with special 
attention directed toward boundary conditions for Jameson’s Multi-grid scheme. 


6.1 Filtered Partial Differential Equations 

The Euler and Navier-Stokes equations may be written in the following general 
form for a non-orthogonal coordinate system: 


_a 

at 



+ 


9F dG 
+ dq 


= s 


( 6 . 1 ) 


Where F and G represent the combined convective and viscous terms, and where S is a 
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source term, which is zero for the unfiltered equations. Consider a general spatial filter 
operator defined as follows: 


I_(f) = f(£,r|) = — 1— f f f(ln) H(^ti) W(f - $,n - n) dfdq (6.2) 

_oo 


where : 


and : 


H(^,r|) = 1 inside the domain 
= 0 outside the domain 

W(^,rj) = Weighting Function 


MM- J 

— oo 


J H(^,Tl)W(^,Tl-Tl)d^dTl 


Apply this filter to Eq. (6.1), 


L 



U ) + iF + 

J * dl ; 


9G 

an 



(6.3) 


(6.4) 


(6.5) 


The first term may be integrated as: 
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1 1 u > W(^-S.n-n-Tl) d^dn (6.6) 

-oo 

oo OO 

J f (r 1 u)H(5,n)W(5-5,n-ii-n)d5dn 
= i|-(Xr T u) 

X dt 


The filtered form of Eq. (6.1) becomes : 


-|-(r 1 U) 


= xJ 

= ii. 

x at 


•j-(J _, U ) + L 

at 

The goal is now to formulate a discrete approximation to this equation which will be 
solved on the first coarse mesh level. Unfortunately, the present equation gives a time- 

dependent equation for J _1 U and not the filtered solution U . It is important to note 
that in general J -l U * J -1 U . Therefore, we will now define the following "average" 
solution, U : 


£ 9G 

a^ an 


= L(S) 


(6.7) 


0 = 


rH j 

F 


In terms of U eq. (6.7) may be written as: 


F 

at 


u) = l 


aF 3G 

as an 


(6.8) 


(6.9) 


In order to put Eq. (6.9) into the same form as the original equation, Eq. (6.1), we define 
filtered approximations to and as follows: 
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= F(U)$ G 11 =G(U) T1 (6.10) 

Then adding these filtered spatial derivatives to each side of Eq. (6.9), the filtered 
equation takes on the same form as the original governing equations (Eq. (6.1)): 

|-(Fu ) + !r + l“r = S (6.1 1) 

at a^ an 

where 


S = L 


S 


aF 


a6 

an 


d¥ BG 

* 


Note that the source term which appears in the filtered equation can not be simplified 
any further since: 


L 

l 


dF(U) 

aG(U) 

an 


,aF(u) 

aG(U) 
* ^ 


(6.12) 


because F and G are non-linear functions of U and the metrics of the coordinate 
transformation, and U * U . 


If we define a discrete approximation to the filter operator (i.e., the fine to coarse 
mesh transfer operator in Multi-grid terminology), Eq. (6.1 1) is the differential equation 
which corresponds to the coarse mesh discrete equation solved in Multigrid 
acceleration schemes To illustrate this point consider the following discrete filter 
operator. Assume a uniform mesh in computational space with = At| = 1 and a 
weighting function, W, defined as having unit influence over the four fine-mesh cells 
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surrounding a given mesh point, and zero influence elsewhere (Figure 6.1): 

W(£ - -ti)= 1 if -1 <£-£< 1 and -1 <r\-r\ < 1 

= 0 if lC-^l > 1 or | t) — Tj | > 1 


n-n 


+i 


4 

3 


+ 

+ 




+ 1 

i + 

+ 2 







Figure 6. 1 . Weighting Function Description 


The corresponding "average" solution, U , is defined as: 


r 1 u _ urUj) 

F Ur 1 ) 

Inside the domain (H = 1 ), with W as defined above: 

ii ii li 

r'u = L(J -1 u) = j J f r'u d^dti = 1 J J u(r‘ dfdri) = 1 J J u dv 


-i -i 


-i -i 


-i -i 


= i£A v iUi 

1=1 


( 6 . 13 ) 


( 6 . 14 ) 


( 6 . 15 ) 


In the last step we have used the mean-value theorem with U; being the cell centered 
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value of the solution on the fine mesh. In a similar fashion : 


l 1 


r> = ur 1 > = 1 f f <r' d^) = | £ av, 

*■ -i -i *■ 

and U is therefore given as: 


(6.16) 


I AVjUi 

u=^ * U 2h =lg h U h (6.17) 

£ av; 

i = 1 

The "average" solution U is exactly the same as the transferred solution, U , used in 
Jameson’s FAS Multi-grid formulation (see Eq. (5.9) ). Similarly: 


dF 3G 

^ an 

= ^r-5— Z(Flux)i 

ZAVi « 

i=l 


-HI 


-i -i 


3F _ 3G 

an 


d^dn 


(6.18) 


= l£ h (F h -L h U h ) 

Which is the same as the transferred fine mesh residual as defined by Jameson (see Eq. 
(5.10)). 

In summary, Eq. (6.11) is the differential counterpart to the discrete coarse mesh 
equations used in Multi-grid acceleration schemes. Expressed in this form, there is a 
direct relationship between the filter operator and the resulting filtered governing 
equations. While the present discussion details a single application of the filtering 
process, if the process is repeated, filtering over increasing length scales, higher and 
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higher levels of the filtered equation can be formulated. Further, each level can be cast 
in the same form as that of eq. (6.11). From each new level of filtered equation the 
corresponding discrete coarse mesh Multigrid equation can be determined. Based on 
this general derivation, it is clear that Jameson’s Multi-grid acceleration scheme is one 
of a large number of Multi-grid formulations which can be constructed from these 
equations. 

6.2 Coarse-Grid Boundary Conditions 

The boundary conditions on the coarse mesh levels are applied at the cell faces, in 
the same manner as the fine-mesh boundary conditions, with the addition of forcing 
terms to allow them to be solved with fine-mesh accuracy. These boundary forcing 
terms will be derived in a similar fashion to the interior forcing terms, using the above 
general Multigrid approach to determine the correct way to transfer the boundary values 
to the coarse-meshes. 

At a boundary, as in Figure 6.2, the function H is zero for the image cells and 1 for 
the interior cells. The weighting function, W, is the same function as for the interior, 


given in Eq. (6.12). 
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Figure 6.2. Coarse-Mesh Boundary Transfer Notation 


Applying the filter, L , to the fine-mesh cells adjacent to the boundary, and integrating 
as before, the transferred boundary value on the coarse mesh, U b , is defined as: 


U b = 


I AV.Uj 

I* l 


(6.19) 


The coarse-mesh boundary conditions are now formulated in terms of the transferred 
solution using the same discrete Multigrid equations as for the interior points: 

L 2h U b = F 2h = Ih h (F h - L h U h ) + L 2h (Ih h U h ) (6-20) 

Since the boundary conditions on each level are satisfied exactly, the quantity, 
l2h(P h _ L h U h ) is always zero, and therefore, the forcing term is : 


F 2hb = L 2h ( ) (6-21) 

which is a function of the transferred boundary value, U b , and the transferred interior 
solution at point 2c (Figure 6.2). As an example, take the inviscid flow tangency 
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condition: 


q Sfc =q^ ; qn b =o (6.22) 

For the tangential velocity, q Sb , the boundary condition on the fine mesh in Multigrid 
notation is: 

LhU h = F h (6.23) 

where 


- q Sb q S2( 

F h =0 

For the coarse mesh levels: 

J-2hU2h = F2h (6.24) 

where 

) 2h 

^ (0) (0) 

F 2h=q s „ -q^ 

and where the quantities in parentheses denote the Runge-Kutta stage, with zero 
representing the initially transferred quantities. The updated boundary value at each 
Runge-Kutta stage is therefore given as: 

(W 00 r (0) (0) , 

q Sb =q s * +[q Sb -q $2c 1 (6.25) 

The conservation variables in the image cells are then updated by linear extrapolation 
from the interior based on the updated value of U b : 

U lc =2U b -0 2c (6.26) 

For the Navier-Stokes equations, the forcing term at a solid wall essentially adds a 
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partial slip condition as the mesh becomes too coarse to fully resolve the boundary 
layer. In the limit, as the mesh becomes very coarse, this is equivalent to applying the 
full inviscid slip condition at the wall. 


CHAPTER 7 - RESULTS 


The numerical algorithm described in the previous sections has been applied to a 
variety of both inviscid and viscous 2-D flow computations in order to determine the 
effect of the Multigrid boundary conditions on convergence. Solutions in each case are 
given to validate the accuracy of the present code, and then convergence histories for 
applications with and without the coarse-mesh boundary forcing terms are shown. It is 
stressed here that all Multigrid and single-grid solutions were the same. References to 
forcing terms refer to those implemented in the boundary cells (Chapter 6) and should 
not be confused with the forcing terms for the interior points given in chapter 5. The 
number of grids specified refers to the number of mesh levels used in the Multigrid 
calculation, with 1 being the finest mesh. 

Convergence is measured in terms of the average residuals of the conservation 
variables which have been normalized at each point by the local time-step. One 
Multigrid Cycle is defined as a complete sequencing through all mesh levels to the 
coarsest mesh and then back to the finest mesh again. For single-grid solutions, this is 
the same as one time-step, or one iteration. 

Computations were performed on a Gould NP-1 mainframe computer, and on an 
Ardent Titan graphics workstation. CPU times on the Titan with a vectorized code 
were equivalent to those on the NP-1 with an unvectorized version. All CPU times are 
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given in minutes. 

7.1 Invisdd Channel Flow 

Initial code validation was performed for subsonic, transonic, and supersonic flow 
in a channel with a circular-arc bump on the lower wall. This is a fairly common 
problem, used in other works as an initial validation tool.^ 2 ^ 13 ^ 15 ^ 27 ^ A 129x33 
algebraically stretched mesh is used for all 3 cases with minimum increments : 
Ax = Ay = .015625. The subsonic and transonic cases use a 10% bump, and the 
supersonic case uses a 4% bump, both of which are shown in Figure 7.1. The channel 
is 3 chords long and 1 chord high. 

7.1.1 Subsonic Case 

Figure 7.2 shows the solutions obtained for the M«„ = .5 case. The spikes in the 
total pressure loss at the leading and trailing edge (Figure 7.2b) are due to the artificial 
dissipation in the numerical scheme activated by the stagnation point pressure 
gradients, and to the discontinuity in the grid metrics there. The trailing edge pressure 
loss is convected downstream and slightly affects the symmetry of the solution. 

Figure 7.3 presents the convergence histories for the subsonic case with and without 
forcing terms for 3, 4, and 5 grid levels. From these plots, it can be seen that the fine- 
grid convergence level of a Multigrid calculation is affected by errors introduced as a 
result of the incorrect formulation of the boundary conditions on the coarse levels. Note 
that with forcing terms applied, all 3 solutions reach the same convergence level as the 
fine-mesh. As the number of grid-levels increases, the coarse-mesh truncation error 
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grows (without boundary forcing terms) due to the increasing mesh coarseness. This is 
summarized in Figure 7.4 for the 3,4, and 5 level solutions without forcing terms. The 
forcing terms therefore permit any number of grids to be used (as permitted by the 
initial mesh resolution of the finest level) while still maintaining fine-mesh accuracy. 

Figure 7.5 shows the efficiency gains realized in the Multigrid calculations for 
various numbers of grids. The single-grid convergence rate is shown and compared to 
2,3,4, and 5-grid rates for Multigrid cycles, and CPU time. CPU time is important in 
Multigrid calculations, since no actual saving are realized by Multigrid acceleration if 
the number of cycles required for convergence, though fewer, take longer to perform. 
Note that the 3,4, and 5 grid solutions use the same amount of CPU time to reach a 
converged solution. This shows that the additional work required above 2 levels is 
negligible for the present Multigrid algorithm, illustrating its efficiency. 

7.1.2 Transonic Case 

Figure 7.6 shows the solutions for the supercritical, M» = .675 case. A small region 
of locally supersonic flow develops over the bump, with a shock occurring over 3 grid 
points at approximately 70% of the chord. The Mach number upstream of the shock is 
1.44, and 5% total pressure loss is generated and convected downstream along 
streamlines. For a normal shock of 1.44, shock tables give a total pressure loss of 
5.24%, which agrees very well with the calculated loss. 

Figure 7.7 shows the convergence histories of the transonic case for solutions with 
and without boundary forcing terms. The same trends exist here as were realized for 
the subsonic case, with the untreated boundary conditions on the coarse-levels 
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introducing coarse-mesh truncation error into the fine-mesh solution. 

Figure 7.8 shows the difference between the single grid and 5-grid convergence 
rates compared to both Multigrid cycles and CPU time. The gains realized by Multigrid 
are somewhat less than for the subsonic case, but this can be attributed to the increased 
convergence rate of the single-grid solution. 

7.1.3 Supersonic Case 

Figure 7.9 presents the results for the M,* = 1.4 channel flow case. Oblique shocks 
form at both the leading and trailing edge and expand as they propagate downstream 
and towards the upper boundary. The leading edge shock angle is such that a standing 
normal shock exists at the upper boundary where the leading edge shock impinges. The 
trailing edge shock is somewhat thicker due to the effect of the expansion waves 
generated over the bump. As the leading edge shock is reflected from the upper wall, it 
is weakened as it passes through the trailing edge shock and then intersects the lower 
wall before being weakly reflected. 

Figure 7.10 shows the single grid convergence rates for the supersonic case 
compared to the 5-grid computations. Applying or not applying the boundary forcing 
terms made no difference in this case. The single grid solution converges very quickly, 
while the Multigrid solution is somewhat slower than in the previous cases, resulting in 
a large decrease in the overall Multigrid efficiency. This is due to the increased flow 
speed. 

As the overall speed of a flow increases, single-grid convergence rates tend to 
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increase, and Multigrid efficiency drops. This can be shown by a comparison of the 
convergence rates, with and without Multigrid acceleration, of the above 3 cases, shown 
in Figure 7.11. The convergence rates of the single-grid solutions increase with Mach 
number, while the Multigrid convergence rates remain approximately the same. 
Therefore, the Multigrid efficiency gains are limited by the fine-mesh. 

7.2 Viscous Flat Plate 

As a first step in validating the present code for the solution of the full Navier- 
Stokes equations, laminar viscous flow was computed over a flat plate. A 65x33 
algebraically stretched mesh was used for these computations, shown in Figure 7.12, 
with minimum increments of Ax min = .01 and Ay min = .0025. The leading edge of the 
plate begins one unit in from the inlet boundary and extends the length of the domain. 

Flow tangency is enforced along the lower wall upstream of the plate. At the upper 
wall, the normal velocity is extrapolated from the interior, and the pressure, tangential 
velocity, and density are specified as having freestream values. The inflow boundary 
conditions are the same as those used in the inviscid channel flow solutions. However, 
due to the non-isentropic nature of the boundary layer that passes through the outflow 
boundary, pressure is specified and the other variables are extrapolated from inside the 
domain. For this calculation, the exit pressure is specified as freestream. A no slip 
condition is enforced along the flat plate, with pressure and density obtained by 
specified adiabatic conditions at the wall, ( 9T/ dt| ) wa n = 0, and zero pressure change 
normal to the wall, ( dP/ 3ti ) W aii = 0- 
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Figure 7.13 shows the solutions for M. = .5 and Reynolds Number: 
"Re x ) x= i = 8000. The skin friction coefficient and boundary layer profiles are compared 
to results obtained from a 2-D compressible boundary layer calculation 1281 utilizing 99 
points normal to the wall. The velocity and temperature are normalized by freestream 

values and plotted verses the boundary layer parameter, tj = y/x^Re x , where Re x is 
based on distance from the leading edge of the plate. 

Figure 7.14 presents the convergence histories for this case with and without 
boundary forcing terms. As in the inviscid computations, the level of error becomes 
higher as the number of grids increases. Note that with the forcing terms applied, the 
same convergence level is reached for both the 4 and 5 grid calculations. For this 
viscous calculation, the error introduced by the coarse-mesh boundary conditions 
without forcing terms seems to be greater than for the previous inviscid cases. This can 
be attributed to the presence of the boundary layer. As the grid levels become coarser, 
boundary layer resolution is not possible and so the no-slip boundary condition on these 
levels is physically incorrect. The forcing terms essentially add a partial slip condition 
as the first mesh line becomes farther away from the boundary, and in the limit of the 
coarsest possible mesh, this is equivalent to specifying flow tangency along the wall. 

However, merely specifying flow tangency on the coarse mesh levels without the 
forcing terms still introduces error into the fine-mesh solution, as can be seen in Figure 
7.15a. When the forcing terms are utilized, all coarse-mesh truncation error is 
effectively eliminated. In Figure 7.15b, this flow tangency condition on the coarse 
mesh levels with forcing terms is compared to the no-slip condition with forcing terms, 
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and the flow tangency convergence is somewhat faster. This may may be due to an 
over-relaxation at the boundaries resulting in increased steady-state convergence. 

Figure 7.16 shows the gains in efficiency realized by Multigrid solutions for the flat 
plate case. They are not as high as those obtained in the previous inviscid cases which 
is typical for viscous solutions. The boundary forcing terms, although reducing the 
error for the 5-grid solution, did not improve the convergence rate as was initially 
hoped. 

7J Viscous Circular Arc Bump 

Subsonic, laminar flow through a channel with a 5% circular-arc bump on the lower 
wall is computed. l ,5 H 271 l 291 a 65x33 mesh is used in these calculations and is shown 
in Figure 7.17. The minimum Ax and Ay are the same as for the flat plate mesh, and the 
stretching is applied at both the leading and the trailing edge. 

Flow tangency is enforced both along the upper wall and the lower wall upstream of 
the leading edge. A no-slip condition is applied over the bump and downstream of the 
trailing edge. This case may also be referred to as a 10% bi-circular arc cascade with a 
sting mounted at the trailing edge and extending downstream. 1151 The normal 
momentum equation (Eq. 4.20) is used to determine the pressure, and density is 
obtained from the wall temperature, specified for this case as the freestream total 
temperature: T wa n = T 0 _ . The inflow and outflow boundary conditions are the same as 
for the flat plate case. 


Figure 7.18 gives the solutions for M«, = .5 and a Reynolds Number (based on 
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chord) of 8000. Boundary layer separation occurs at 81% of the chord, with the 
separation bubble extending downstream to 120% chord, where the flow re-attaches. 
Skin friction and static pressure results agree very well with those given by 
KallinderisJ 29 ' The velocity profiles agree fairly well with those reported by Rhie.^ 
However, a slightly thicker boundary layer is predicted by Rhie, possibly due to higher 
levels of artificial dissipation, and the discrepancy increases with distance over the 
bump. The separation bubble occurring at the trailing edge may be seen in the contour 
plots and velocity vector plots. 

Figure 7.19 shows the Multigrid results for this case with and without boundary 
forcing terms for 3, 4, and 5 grids. The increase in error with number of grids is greater 
than that for the flat plate case, possibly due to more complex flow phenomena and the 
additional stretching at the trailing edge. 

Figure 7.20 summarizes the increased effect of the coarse-mesh truncation error 
(without forcing terms) on the fine-mesh as the number of mesh levels increases. The 
effect of applying flow tangency on the coarse grids for this viscous case is shown, with 
and without forcing terms, in Figure 7.21a. The results are similar to the corresponding 
flat plate case (Figure 7.15a), although there is no apparent increase in convergence rate 
for the flow tangency condition with forcing terms (Figure 7.21b). 

Figure 7.22 presents the efficiency of the Multigrid calculations compared to both 
Multigrid cycles and CPU time with the single-grid convergence. The increase in 
efficiency measured by CPU time is less than that for Multigrid cycles due to the 
additional work required in calculating the viscous stresses. 
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7.4 VKI Gas Turbine Rotor Blade 

In the interest of applying the present algorithm to a more realistic problem, inviscid 
and viscous calculations are performed here for a VKI gas turbine rotor blade. 1301 The 
common geometric parameters for cascade calculations are defined in Figure 7.23, and 
for this case are specified as: 

— = .697 y= 33.35 P 1 =24° 

c 

Where g/c is the gap-to-chord ratio, y is the blade stagger angle, and is the inlet flow 
angle. 

The present computations were performed for M! = .19 and M 2i< = 1.0, where M^ 
is the specified outlet isentropic Mach number which sets the ratio of exit static pressure 
to upstream total pressure (see Eq. 4. 19) This ratio is .53 for the current value of . 

The flow is turned 96 degrees to a final steady-state outlet flow angle of approximately 
72 degrees. 

7.4.1 Inviscid Case 

Figure 7.24 shows the 73x17 C-mesh used for the inviscid calculations. At the 
inlet, Riemann Invariant boundary conditions are used (see chapter 4) with specified 
total pressure, total temperature, and flow angle. The outflow boundary conditions use 
characteristic variables with a specified exit pressure. Periodicity is enforced along the 
upper and lower boundaries, and flow tangency is enforced on the blade. 

Figure 7.25 gives the inviscid results. All surface quantities are plotted along the 
blade’s coordinate system, not the axial distance. The local isentropic Mach number is 
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given as : 



2 

1 - 1 


Mi s — 

( Y- 1 ) 

(p/Po,) <Y -' Vr 



Comparisons are made with blade-to-blade H-mesh calculations performed by Arts' 
which used a corrected viscosity scheme, and with experimental results obtained by 
VKI. 

The flow is accelerated on the suction surface until a weak shock forms. The 
present calculations have smeared this region, and have underpredicted the flow 
expansion, possibly due to the coarseness of the mesh and the type of artificial 
dissipation that is used. The leading edge pressure gradients are very high, activating 
the 2nd order dissipative terms and creating large total pressure losses. Large total 
pressure losses are also incurred at the trailing edge due to the wedge-type truncation 
verses a smooth, blunt rounding of the trailing edge. In order to realize minimum total 
pressure loss for inviscid cascade computations, it is essential that the normal 
momentum equation, tangential momentum equation, and flow tangency all be satisfied 
simultaneously at each point along the blade. 191 However, this method was not used in 
the present work. The weak suction surface shock and supersonic bubble may be seen 
in the Mach contours at about 50% chord. 

Figure 7.26 shows the convergence histories for the inviscid case. Since the 
coarse-mesh truncation error is quite low for only 3 levels, the boundary forcing terms 
had no effect on the convergence levels for this case. Due to a large portion of high- 
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speed flow, the single-grid convergence was quite rapid, and Multigrid therefore did not 
produce large savings. 

7,4.2 Viscous Case 

Figure 7.27 shows the 129x41 C-mesh that is used in the viscous computations. 
Results for laminar flow at a Reynolds Number of 10,000 are presented in Figure 7.28. 
The isentropic Mach number compares quite well with Arts’ f301 computations and the 
experimental results J 30 * However, the shock on the suction surface remains smeared 
and the flow upstream of it under-expanded due to viscous smearing which occurs in 
regions of high adverse pressure gradients. A large portion of low-speed flow exists 
near the stagnation point, seen by the delayed onset of boundary layer development on 
the pressure surface (Contour Plot). Laminar separation can be seen occurring at about 
40% of the blade chord on the suction surface (velocity vectors). 

Figure 7.29 shows the convergence histories for this case. 1 and 3 grid convergence 
rates are plotted versus Multigrid cycles and CPU time. Multigrid produces substantial 
savings for this case. 




Figure 7.1. 129x33 Computational Grids For Inviscid Channel Flow, (a) 10% Circular 
Arc Bump (b) 4% Circular Arc Bump 
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Figure 7.2. Mo. = .5 Inviscid Channel Flow Solution, (a) Surface Mach Number (b) 
Surface Total Pressure Loss (c) Mach Contours, AM = .05 (d) Total 
Pressure Loss Contours, AP 0 = .01 
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Figure 7.2. Continued. 
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Figure 7.3. Effect of Boundary Forcing Terms on Inviscid M** = .5 Convergence. 3 
Grids: (a) No Forcing Terms (b) Forcing Terms (c) Summary. 4 Grids: (d) 
No Forcing Terms (e) Forcing Terms (0 Summary. 5 Grids : (g) No 
Forcing Terms (h) Forcing Terms (i) Summary 
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Figure 7.3. Continued. 
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Figure 7.3. Continued. 
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Figure 7.4. Increase In The Influence Of Coarse-Mesh Truncation Error On Fine-Mesh 
Convergence As The Number Of Grids Increases 
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Figure 7.5. Multigrid Efficiency For Inviscid, = .5 Case, (a) 1 Grid (b) 2-5 Grids 
(c) 1-5 Grids Multigrid Cycle Comparison (d) 1-5 Grids CPU Time 
Comparison 
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Figure 7.6. = .675 Inviscid Channel Flow Solution, (a) Surface Mach Number (b) 

Surface Total Pressure Loss (c) Mach Contours, AM = .05 (d) Total 
Pressure Loss Contours, AP 0 = .01 
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Figure 7.6. Continued. 
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Figure 7.7. Effect of Boundary Forcing Terms on Inviscid M,. = .675 Convergence. 4 
Grids: (a) No Forcing Terms (b) Forcing Terms (c) Summary. 5 Grids: (d) 
No Forcing Terms (e) Forcing Terms (f) Summary 
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Figure 7.7. Continued. 
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Figure 7.7. Continued. 
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Figure 7.8. Multigrid Efficiency For Inviscid, M.. = .675 Case, (a) 1 Grid (b) 5 Grids 
(c) 1 and 5 Grids Multigrid Cycle Comparison (d) 1 and 5 Grids CPU Time 
Comparison 
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Figure 7.9. M» = 1.4 Inviscid Channel Flow Solution, (a) Surface Mach Number (b) 
Surface Total Pressure Loss (c) Mach Contours, AM = .053 (d) Total 
Pressure Loss Contours, AP„ = .0052 
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Figure 7.9. Continued. 
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Figure 7.10. Multigrid Efficiency For Inviscid, M.. = 1.4 Case, (a) 
(c) 1 and 5 Grids Multigrid Cycle Comparison (d) 1 
Time Comparison 
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Figure 7.1 1. Effect Of Mach Number On Convergence For Inviscid Channel Flow, (a) 
1 Grid (b) 5 Grids 
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Figure 7.12. 65x33 Flat Plate Computational Mesh 



Figure 7.13. Viscous Flat Plate Solution; M*. = .5, Re x ) x=1 = 8000. (a) Skin Friction 
Coefficient (b) Velocity Vectors (c) Velocity and Temperature Profiles at 
x=. 25 (0 Velocity and Temperature Profiles at X=1.25 
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Figure 7.14. Effect of Boundary Forcing Terms on Viscous Rat Plate Convergence. 4 
Grids: (a) No Forcing Terms (b) Forcing Terms (c) Summary. 5 Grids: (d) 
No Forcing Terms (e) Forcing Terms (f) Summary 
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Figure 7.14. Continued. 
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Figure 7.15. Flat Plate Convergence With Flow Tangency Specified On The Coarse 
Levels, (a) With and Without Forcing Teims (b) Flow Tangency With 
Forcing Terms versus No-Slip Boundaries With Forcing Terms 
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Figure 7.16. Multigrid Efficiency For Viscous Flat Plate, (a) 1 Grid (b) 5 Grids (c) 1 
and 5 Grids Multigrid Cycle Comparison (d) 1 and 5 Grids CPU Time 
Comparison 







Figure 7.17. 65x33 Computational Mesh for Viscous Flow Over a 5 % Circular Arc 
Bump 
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Figure 7.18. Viscous Channel Flow Solutions Over a 5 % Circular Arc Bump, (a) Skin 
Friction Coefficient (b) Normalized Static Pressure Distribution. Velocity 
Profiles: (c) X=1.25 (d) X=1.50 (e) X=1.75. (f) Mach Contours, AM = .05 
(g) Total Pressure Loss Contours (h) Trailing Edge Velocity Vectors 











Figure 7.18. Continued. 
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Figure 7. 19. Effect of Boundary Forcing Terms on Viscous Circular Arc Convergence. 

3 Grids: (a) No Forcing Terms (b) Forcing Terms (c) Summary. 4 Grids: 
(d) No Forcing Terms (e) Forcing Terms (f) Summary. 5 Grids: (g) No 
Forcing Terms (h) Forcing Terms (i) Summary 
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Figure 7.19. Continued. 
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Figure 7.19. Continued. 
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Figure 7.20. Increase In The Influence Of Coarse-Mesh Truncation Error On Fine- 
Mesh Convergence As The Number Of Grids Increases 
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Figure 7.21. Viscous Bump Convergence With How Tangency Specified On The 
Coarse Levels, (a) With and Without Forcing Terms (b) How Tangency 
With Forcing Terms versus No-Slip Boundaries With Forcing Terms 
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Figure 7.22. Multigrid Efficiency For Viscous Circular Arc. (a) 1 Grid (b) 1, 4 and 5 
Grid Multigrid Cycle Comparison (c) 1, 4 and 5 Grid CPU Time 
Comparison 
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Figure 7.24. 73x17 Computational Mesh for Inviscid Turbine Calculations, (a) Global 
View (b) Leading Edge Detail (c) Trailing Edge Detail 
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Figure 7.25. Inviscid Turbine Solutions, (a) Blade Isentropic Mach Number (b) Blade 
Total Pressure Loss (c) Coefficient of Pressure (d) Mach Contours, 
AM = .04 (e) Inflow Velocity Vectors (f) Trailing Edge Velocity Vectors 
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Figure 7.25. Continued 
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Figure 7.26. Convergence Histories for Inviscid Turbine (a) l Grid (b) 3 Grids (c) 1 
and 3 Grids Multigrid Cycle Comparison (d) 1 and 3 Grids CPU Time 
Comparison 
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Figure 7.27. 129x41 Computational Mesh for Viscous Turbine Calculations, (a) Global 
View (b) Leading Edge Detail (c) Trailing Edge Detail 
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Figure 7.28. Viscous Turbine Solutions. Re L = 10.000. (a) Blade Isentrop.c Mach 
Number (b) Blade Coefficient of Pressure (c) Mach Contours, AM - W 
(d) Total Pressure Loss Contours (e) Leading Edge Velocity Vectors (0 
Separation Point Velocity Vectors (g) Trailing Edge Velocity Vectors 
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Figure 7.28. Continued. 
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Figure 7.28. Continued. 
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Figure 7.28. Continued. 
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Figure 7.29. Convergence Histories for Viscous Turbine (a) 1 Grid (b) 3 Grids (c) 1 
and 3 Grids Multigrid Cycle Comparison (d) 1 and 3 Grids CPU Time 
Comparison 
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CHAPTER 8 - CONCLUSIONS 

8.1 Major Contributions and Summary of Results 

This work has presented initial research into the effects of coarse-mesh boundary 
conditions on the convergence of Multigrid Acceleration. An explicit Multigrid 
algorithm has been written and validated for a variety of both inviscid and viscous flow 
computations. The basic algorithm has been fashioned after Jameson’s finite-volume 
Multigrid scheme which utilizes explicit, Runge-Kutta time-stepping. Forcing terms ; 
have been derived and added to the coarse-mesh boundary conditions which permits 
them to be solved with fine-mesh accuracy, i.e., without the coarse-mesh truncation 
error polluting the fine mesh solution. 

In order to derive the correct interpolation procedure, or transfer operator, for the 
solution at boundary points, a new, general approach to formulating the governing 
equations on the coarse mesh levels has also been presented. In this approach, the 
equations on the coarse-meshes have been viewed as a filtered sub-set of the fine-mesh 
equations, and a general filtering operator has been derived which models this filtering 
process. This approach is based on the fact that certain information is not resolved on 
the coarse-mesh due to the increased mesh spacing, and, therefore, is "filtered out". By 
applying the filter to the governing fine-mesh equations, a formal description between 
these and the coarse-mesh equations is given. Then, by specifying a discrete form of 
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the filter, the discrete coarse-mesh equations may be formulated. The advantage to this 
formulation is that any number of Multigrid schemes may be formed by specification of 
different discrete filters. In this work, the filter was chosen such that Jameson’s 
Multigrid scheme was obtained. Then, the boundary transfer operator and coarse-mesh 
boundary forcing terms were derived. 

In summary, the major contributions that have been presented in this work are: 

— A new, general approach to obtaining the coarse-mesh governing equations in 
partial differential form, where they are formally related to the fine- mesh 
equations 

— Formulation of a boundary transfer operator that is consistent with the interior 
scheme 

— Derivation of forcing terms for the coarse-mesh boundary conditions 

— Demonstration of the ability of these forcing terms to allow the coarse-mesh - 
boundary conditions to be applied on any number of mesh levels with fine-mesh * 
accuracy 

Flow calculations were performed for inviscid channel flow over a circular-arc 
bump for subsonic, transonic, and supersonic speeds, a viscous flat plate, viscous 
subsonic channel flow over a circular- arc bump, and inviscid and viscous flow over a 
VKI gas turbine rotor blade. Good agreement with other results was seen in all cases, 
with only a few minor discrepancies occurring. The discrepancies which occurred in 
the viscous bump velocity profiles (Figure 7.18c-e) and the inviscid turbine isentropic 
Mach number distribution (Figure 7.25a) seemed to be related to the artificial 
dissipation model. 

Implementing the coarse-mesh boundary conditions without forcing terms 
introduced coarse-mesh truncation errors into the fine-mesh solutions. This was 
apparent in all cases except the supersonic bump and the VKI turbine blade. For these 
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2 cases, the truncation- error was apparently lower than the local round-off error. The 
error introduced into the fine-mesh solution increased as the number of mesh levels 
increased (Figures 7.4,7.20) due to the larger mesh spacings on the coarsest mesh. 
With forcing terms, however, the fine-mesh accuracy was retained for any number of 
mesh levels used. This becomes very important as problems increase in size and more 
mesh-levels are used. The problems here have used relatively coarse initial meshes 
with few mesh levels in the Multigrid cycle. In larger applications which don’t use 
forcing terms in the boundary conditions, the number of mesh levels that could possibly 
be used would be restricted due to the error that is introduced into the fine mesh 
solution. This would consequently reduce the efficiency of the Multigrid solution. . 
Therefore, the forcing terms, by permitting the maximum number of mesh levels to be - 
used, increase Multigrid efficiency even though they don’t actually increase 
convergence rates. 

The forcing terms along the solid boundaries for the viscous problems were stated 
(Chapter 6) as applying a partial slip condition as the mesh becomes too coarse to 
accurately resolve the boundary layer. As the mesh spacing scale on the coarsest mesh 
level exceeds the length scale of the boundary layer, this should be equivalent to 
enforcing flow tangency along the surface. However, when flow tangency is explicitly 
applied, it was shown that the forcing terms are still required to produce fine- mesh 
accuracy (Figures 7.15, 7.21). 

The efficiency of Multigrid acceleration was seen to decrease with increasing Mach 
number for the inviscid channel flow cases. However, this may be attributed to the 
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increase in efficiency of the single-mesh solutions as Mach Number increases (Figure 
7.11a). The Multigrid convergence rates for the 3 cases were approximately the same, 
with slightly reduced performance in the supersonic case (Figure 7.11b). The efficiency 
of Multigrid acceleration for the flat plate, viscous bump, and viscous turbine blade 
problems was greater in terms of Multigrid cycles than in terms of the actual work 
performed, measured in CPU time (Figures 7.16c,d 7.22c, d 7.29c,d). This is typical of 
viscous problems, and is due to the additional work required to evaluate the viscous 
fluxes. 

8.2 Future Work 

The general Multigrid approach that was formulated may be used to derive any " 
number of discrete Multigrid schemes based on a specified filter, or prescribed method 
of transferring the solution and residuals. In this work, only one such filter was looked 
at: the simple 4-cell volume weighted averaging as given in Jameson’s scheme. Future 
work may entail the derivation of various other filters and a study of their effect on the 
performance of Multigrid acceleration. 

Further study could also be performed into the reasons behind the effect of Mach 
number on convergence, based on the observations made for the inviscid channel flow 
problem. For the Multigrid solutions to this case, not much difference existed between 
the convergence rates, although the supersonic Multigrid convergence was the slowest 
of the three, exactly opposite from the single-mesh convergence rates. 

Finally, a more robust dissipation model could be pursued. Some difficulties are 
realized when performing Euler calculations on very fine, highly stretched meshes. 
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possibly due to the scaling that is used. Alternative scaling theories could be studied, 
i.e., based on something other that the maximum spectral radii of the convective 
Jacobian matrices. 
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